Continuous Time Random Walks: Simulation of continuous trajectories 
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Continuous time random walks have been developed as a straightforward generalisation of classical 
random walk processes. Some 10 years ago, Fogedby introduced a continuous representation of these 
processes by means of a set of Langevin equations [H. C. Fogedby, Phys. Rev. E 50 (1994)]. The 
present work is devoted to a detailed discussion of Fogedby's model and presents its application for 
the robust numerical generation of sample paths of continuous time random walk processes. 

PACS numbers: 05.40.Fb,45.10.Hj, 02.60.Cb 



I. INTRODUCTION 

The analysis of stochastic processes by Bachelier, Ein- 
stein and Langevin [l], [2, Q extensively has inspired sci- 
entist in the last century. First, Fokker-Planck equations 
describing the time evolution of probability density func- 
tions (pdfs) emerged. In addition, mathematical meth- 
ods for proper interpretation of the associated Langevin 
equations have been developed, that allow access to the 
trajectories of individual particles. An intrinsic feature 
of these processes is, that they obey Markov properties 
Therefore, two point statistics are sufficient for a 
complete description of these processes. 

In recent years, processes exhibiting anomalous diffu- 
sion, {x'^{t)) ~ with ^ ^ 1, increasingly have attracted 
attention [5|. Such processes typically are realised in 
complex environments such as porous and disordered me- 
dia, see e.g. Iff] and references therein. In contrast to or- 
dinary diffusion, Markov properties do not hold for these 
processes. Therefore, multipoint joint statistics have to 
be considered for proper description of the dynamics. 
Likewise, two alternative approaches to these processes 
have been evolved. On the one hand, these processes can 
be described by means of fractional Fokker-Planck equa- 
tions, that contain fractional derivatives with non-local 
character. On the other hand. Continuous Time Ran- 
dom Walk (CTRW) processes [3] have been proposed for 
the analysis of the microscopic properties of anomalous 
diffusion processes [8|]. In general, they are specified by 
the iterative discrete equations 5, 9] 
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where {r]i,Ti) is a set of random numbers drawn from 
the pdf ^(?7, r), that vanishes for negative values of r for 
reasons of causality. Frequently, equations ([T]) are used 
to model time-continuous processes with the additional 
assignment 0,Q 



with ti < t < t. 
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With the aid of CTRWs, limiting behaviour and ensemble 
statistics of anomalous diffusion processes become acces- 
sible through Monte-Carlo simulations. For details, the 



reader is referred to recent works by Dentz et al. [g], 
Heinsalu et al. [l^l and Gorenflo et al. [Tl| . 

Recently, Fogedby formulated a continuous description 
of CTRWs, that is based of a set of stochastic differential 
equations p^ . 
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Here, the discrete variable i of equations (H]) is generalised 
to the continuous variable s, that can be associated with 
an intrinsic time of the CTRW. A number of publications 
addressed statistical properties of trajectories of this ap- 
proach [l^ [ij, [l^. It is the aim of the present work, 
to present a robust algorithm for the generation of con- 
tinuous sample paths from Fogedby's equations. By this 
means trajectories can be obtained, that properly exhibit 
the anomalous dynamics of CTRWs on any time scale. 

This work is structured as follows. In the next section, 
some general remarks are made on the definition of con- 
tinuous CTRWs, equations Section |llT] is dedicated 
to the properties of the process t{s), equation pb)) . that 
typically is driven by Levy noise. The numerical sim- 
ulation of continuous CTRW processes is presented in 
section IIVI and exemplified by means of some results in 
section|Vl We conclude with section lVll that summarizes 
our results and suggests future applications. 



II. CTRWS IN THE SPIRIT OF FOGEDBY: 
SOME REMARKS 

First of all, the character of the distributions of the 
random variables 77 for the jump length and r for the 
waiting time has to be addressed. In case of the discrete 
definition, equations ([T]) , a broad class of distributions is 
feasible for this purpose. The continuous Langevin for- 
mulation of Fogedby, equations ([3]), however, requires the 
associated distributions to be stable in order to be prop- 
erly defined. For a short introduction into the concept 
of stable distributions we refer to [5]. In the long time 
limit, however, both approaches are equivalent. 

In 1994, Fogedby introduced the stochastic differential 
equations ([3]) for CTRWs as the continuum limit of the 
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FIG. 1: Examples for the fully skewed Levy distribution ac- 
cording to equation (O for different characteristic exponents 
a. For a — > 1, the pdf converges to S{x — 1). 



path parameter or arc length s along the trajectory [12] ■ 
He considered independent random variables 77(5) and 
t(s) with power law behaviour, that is 



*(r/,r)--ry ^ ^"t ^ for ?7,t>1 



(4a) 



For reasons of normability, Fogedby used cutoffs at low 
values of \ri\ and t, respectively. He mainly was interested 
in the properties of the process x{t), that can be obtained 
from equations ^ by inversion of the latter process. By 
means of his modified power-law distributions, the long- 
time behaviour of processes with power law waiting jump 
length and waiting time distributions could be derived. 
In this context, first properties of the inverse process s{t) 
of equation ([3b| have been addressed. 

Baule et al. investigated the properties of the inverse 
process s{t) in greater detail |13j . In particular, multi- 
time joint probabilities could be calculated. The waiting 
times T were considered to obey one-sided Levy distribu- 
tions with tail parameter a, that are discussed in greater 
detail in the following section. In absence of external 
force terms (F = in equation (|3ap ) analytical expres- 
sions for correlations functions could be derived by ap- 
plication of the inverse Fourier- and Laplace-transforms. 
Recently, these results could be extended to Ornstein- 
Uhlenbeck like processes with a linear repelling term, 
thus F{x) = --fx [H. 

Both works, however, focus on the ensemble statistics, 
whereas the Langevin approach of Fogedby, that is inter- 
esting itself, has not attracted much attention yet. 



III. PROPERTIES OF t{s) AND ITS INVERSE 
PROCESS s{t) 

In this section, we concentrate on the process t{s). Due 
to causality, this process has to be strictly monotonically 
increasing. The increments dt{s) = t{s-\-ds) — t{s), there- 
fore, have to obey fully skewed stable distributions. 



Generally, stable distribution are assigned to the a- 
stable probability distribution functions (pdfs). a-stable 
Levy distributions typically are not available in a closed 
form but are only accessible by means of their Fourier 
transform. We especially consider the distribution 



(5) 

Here, i is the imaginary unit and < a < 1 the stability 
index, that specifies the asymptotic behaviour La{x) ^ 
^.-(i+q) at a: 3> 1. This pdf complies with a common 
parametrization of a-stable pdfs [H, [l3| , 
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.2 a7r^(-l/2) 



for the parameters /3 = 1 and c = (l + tan^ 
Some examples for this distribution are depicted in figure 
[T] An important feature of this representation is, that it 
according to equation ([5]) is defined even for a = 1 by 
means of La=i{x) = 6{x — 1). From the theory of stable 
processes it follows, that the increment dt has to obey 
the distribution ds^°-La{dt/ds°-) [l^. The case a = 1 
with dt = ds complies with a stochastic processes, that 
solely can be described by equation (|3ap with s = t. 

An intrinsic feature of Levy processes is, that trajecto- 
ries contain finite jumps in terms of discontinuities with a 
probability greater than zero. They are continuous only 
from the right ^8'], that is 



lim t(s + As) = t(s) 



(7) 



A sample of a fully skewed Levy process with character- 
istic exponent a = 0.9 is depicted in the upper panel of 
figure m Due to these jumps, the range of values t{s) 
does not cover the full interval [0, 00 [. Rather, the range 
can be specified by means of a set of intervals. Due to the 
monotonic increase of the process t{s) a inverse process 
s{t) exists. This process has been applied for construc- 
tion of the process x{t) = x{s{t)) in the past. It, however, 
a priori is not properly defined for t G [0,oo[ due to the 
jump characteristic of the Levy process. A sample of the 
inverse process is depicted in the lower panel of figure [51 
In general, the meaning of the inverse function has to 
be specified in order to be properly defined on t e [0, 00 [. 
App ropriate definitions for the inverse function e.g. are 

m 

s{t) := inf {s : i(s) > i} or (8a) 
s{t) := sup{s : t{s) < t} . (8b) 

If the process x{s), is steady, these definitions are 
equivalent in the limit As — > for at t € [0,oo[: Since 
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FIG. 2; Trajectories of the process t{s) (upper panel) and 
the associated inverse process s{t) (lower panel) simulated 
for a = 0.9. The finite jumps of the process t{s), that are 
characteristic for Levy processes, are evident. Due to these 
jumps the inverse process s{t), however, a priori is not defined 
on [0, oo[. 



x{s) is steady, 



lim x(s — As) — lim x(s + As) — x(s) (9) 



is valid by definition. Due to the monotonia character 
of the process t{s), t{s) < t{s + As) for As > 0. Con- 
sequently, s < s{t') < s + As if t(s) < t' < t{s + As). 
Then, 



lim x(s) 



x{s{t')) = x{s + As) 



(10) 



is valid for i(s) < t' < t{s + As). Here, for the process 
t{s) only the continuity from the rights has been used. 
For steady processes x{s) in the limit of infinitesimal As, 
thus, no additional definition for proper interpretation of 
the inverse function is required. 

If unsteady jumps of x{s) coincide with those of t(s), 
the latter argumentation is not feasible. We, however, 
restrict to stochastic processes x{s) with Gaussian noise, 
that are steady with probability 1 for s € [0, oo[. 



IV. NUMERICAL SIMULATION OF SAMPLE 
PATHS 

An equivalent formulation of equations ^ is given by 
the integral equations 



x{s) = x{0) + J ds'F{x{s')) + J dWis') (11a) 





t{s) = i(0) + J dL^{s') 





(lib) 



where dW and dL^ are the infinitesimal increments of 
Wiener and a-stable Levy processes, respectively. For 
numerical integration these equations have to be discreti- 
sized with an adequate discrete increment As. Appli- 
cation of the Euler scheme for numerical evaluation of 
equations pTjl then yields 

a;(s-|-As) = x{s) + AsF{x{s)) + ri{s, As) (12a) 
t{s + As) = t{s) + Tais,As) . (12b) 

Here, the random variables r]{si, As) independently have 
to be drawn from a Gaussian pdf with variance — As. 
The variables Ta{si, As) have to comply with the distri- 
bution -^^La{-^js-)- The efficient numerical generation 
of these random numbers is addressed in appendix [A] 
Due to the absence of forcing and the purely additive 
character of the noise, the integration of t(s) by means 
of the Eulcr scheme is exact. For numerical integration 
of the process x(s) with Gaussian noise, alternatively ad- 
vanced discretization schemes can be applied ,19J . 

For numerical simulation of trajectories x{t) at discrete 
times tj := j At, j = 0,. . . ,N, the inverse s{t) does not 
have to be calculated explicitly. Instead, the following al- 
gorithm can be applied, that incorporates definition (j8ap 
for the inverse process: 

• Initialisation of Xs{0) and ts{Q), set s = 

• for every j = to N: 

1. while {ts{s) <tj): 

(a) calculate Xs{s + As) and ts{s + As) from 
eqns. (fT2| 

(b) increase s by As 

2. set x{tj) :— Xs{s) 

The discretization of t, At, is given by the desired sam- 
pling rate of the simulated process. The optimal value for 
the discretization of the intrinsic variable s. As, depends 
on characteristic length scales of the process x{s) and the 
desired accuracy of the resulting process x{t). Typically, 
As has to be adjusted, such that the right hand side term 
of the discrete equation, Asi^(a:(si)) -I- 77(si, As), with a 
sufficient probability is less than the desired accuracy of 
the simulation. This argument is illustrated within the 
scope of the following section. Too small values for As, in 
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FIG. 3: Determination of the intrinsic increment As. From 
the transition pdf of Ornstein-Uhlenbeck processes, that is 
available in a closed form, the square root of the means square 
deviation of the increment Aa; as function of the increment As 
can be derived. From inspection of this graph, As = 0.0001 
seems to be sufficient for the current purpose. 

turn, may bias the accuracy of the numerical evaluation 
of equations due to discretisation errors and reduce 
the speed of the algorithms. The choice for As therefore 
is a trade-off between accuracy of the discretization, va- 
lidity of the inversion of the process t{s), computer time 
and discretization errors. For details concerning the nu- 
merical evaluation of the discretized equations (fT2|) . the 
reader is referred to the book of Kloeden and Platen [l^ . 

V. EXAMPLES 

For exemplification of the simulation procedure and 
characteristic properties of continuous CTRWs, processes 
with F{x) — —X are considered. The process x{s) then 
is an ordinary Ornstein-Uhlenbeck process 

dx{s) = -jxdt + v^dW{s) (13) 

with £) = 7 = 1. For later comparison with analytical 
results, a;(0) = 1 is used as starting value for the process 

For Ornstein-Uhlenbeck processes, joint pdfs for finite 
time increment can be calculated in a closed form [J. 
Then, the statistics of the increments Ax := x{s + As) — 
x{s) can be considered as a function of the discretization 
As. The distribution of Aa; as a function As is Gaussian 
with variance 

((Aa;)2)^2(l-e-^^) . (14) 

In order to estimate an appropriate value for As, the root 
mean square deviation has been investigated. From the 
inspection of figure [3l As — 0.0001 has been selected for 
application in the numerical procedure. The maximum 
deviation of the definitions ([5]) from one another is ~ 
10~^, which is accurate enough for the current purpose. 



10000 data points with time increment At = 0.001 
have been generated for several values of a. The trajec- 
tories of the respective processes are exhibited in figure 
m From the sample paths, the influence of the subor- 
dinating process t{s) on the dynamics becomes evident. 
For a = 1 an Ornstein-Uhlenbeck process is recovered. 
With decreasing a waiting events start to dominate the 
process. This also can be seen from the evolution of the 
increment pdfs, that is depicted in figure [5] 

Recently, the fractional extension of Ornstein- 
Uhlenbeck processes has been investigated by Baule et 
al. ^5| , starting from the fractional Fokker-Planck equa- 
tion for the time evolution of ensembles of particles. For 
the Ornstein-Uhlenbeck process x(s) with initial value 
a::(0) = 1, that has been considered in this section, for 
^2 > ti eventually an analytical expression for the corre- 
lation function could be derived, 

4-a ^ ( 4-a.\n 

{x{h)x{h)) = y: ^ 'i^, (15) 

X2Fi(^a,-an,a + l-,^-^^ +Ea{-t^) . 

Here, F denotes the Gamma function, the one- 
parameter Mittag-Lefiler function and 2^1 (a, b, c; z) the 
Gaussian hypergeometric function. For details see the 
references provided by Baule et al. p^ . 

On the other hand, the correlations can be estimated 
from ensemble averages of simulated trajectories x{t). A 
comparison of the analytical result by Baule et al. with 
the correlations estimated from our simulated trajecto- 
ries is exhibited in figure [6l Perfect coincidence of these 
two approaches is observed, that never could be com- 
pared before. 

VI. SUMMARY AND CONCLUSIONS 

A method for the accurate and efficient simulation 
of continuous trajectories of Continuous Time Random 
Walk (CTRW) processes has been proposed, that relies 
on the representation through Langevin equations pro- 
posed by Fogedby. It is based on the simultaneous simu- 
lation of two stochastic processes, one of which is driven 
by Levy noise. 

Within the scope of section Hill the construction of the 
process x{t) = x{s{t)) with the aid of the inverse s{t) of 
the Levy process t{s) has been discussed in great detail. 
The unique existence of the inverse of t(s) for any t typ- 
ically has been assumed in the past 13]. However, the 
meaning of the inverse process in fact has to be specified 
in detail at discontinuities of t{s) in order to guarantee for 
unique existence, see e.g. equations ([8|). We would fike to 
emphasize, that these additional specifications infiuence 
neither the trajectories x{t) nor the ensemble statistics, 
if the trajectories x{s) are continuous with probability 1. 
In the subsequent sections, we mainly focussed on this 
specific case. 
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a=1.00 



a=0.85 




a=0.95 




a=0.90 





a=0.80 




a=0.50 




FIG. 4: Sample trajectories of CTRWs with linear repelling force F{x) = —x for different stability indices a. The solid line 
corresponds to the process x{t) whereas the dashed line indicates the corresponding s{t). a = 1 complies with the ordinary 
Ornstein-Uhlenbeck process. With decreasing a the process is dominated by waiting events indicated by constant x and s. 



Comparison with recent analytical results by Baiile et 
al. [iBl has been used for validation of the simulation 
procedure and showed compliance of the results. Due 
to the non-Markovian character of fractional processes, 
higher order joint statistics are of great interest [20]. We 
propose the use of probabilistic methods for numerical 
calculation of these functions. 
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FIG. 5: Increments distributions for Aa; := x{t+T) — x{t) with 
T = 0.001 for some processes depicted in figure|4l For reasons 
of clearness the individual pdfs are shifted in vertical direc- 
tion by a constant factor. It is evident, that with decreasing 
stability index a a central peaks evolves, that corresponds to 
persistent regions due to waiting events. On the other hand, 
the distributions broaden indicating a higher probability of 
the occurence of extreme increments. 
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FIG. 7: Performance of the generation of skewed a-stable 
random numbers. The solid line exhibits the pdf obtained 
by numerical integration of ([5} for a = 0.8. For a; 3> 1 the 
pdf shows power law decay with the exponent 1 -I- 0.8, that 
is depicted dashed. The points mark the pdf obtained from 
a sample of 10^ random numbers, that have been generated 
by means of equation (|A3|I for the same stability index. The 
analytical pdf evidently is well reproduced by the sample of 
random numbers. 
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FIG. 6: Correlation function of the process depicted in fig- 
ure |4] for a — 0.8. The solid line marks the analytical so- 
lution (|15p derived by Baule et al. [isl ]. The circles indicate 
the correlation obtained from the analysis of an ensemble of 
500000 trajectories. Since a perfect coincidence is observed, 
this evaluation is proposed as a benchmark for accuracy of 
the numerical implementation. 



APPENDIX A: GENERATION OF RANDOM 
VARIABLES WITH SKEWED a-STABLE PDF 
ACCORDING TO EQUATION ^ 

Skewed Levy-stable random numbers efficiently can be 
generated by means of the algorithm proposed in 
We adapted this algorithm to our definition of the skewed 
Levy distributions, ©■ 

The random numbers Ta{si^As), that are required 
for numerical integration of the process t{s), then for 
< a < 1 efficiently can be generated by means of this 



• Generate a random variable Vi uniformly dis- 
tributed on ] — Tr/2, 7r/2[ and an independent expo- 
nential random variable Wi with mean 1. Several 
optimized random number generators are available 
for this purpose. In case of doubt, V and W can be 
obtained from two independent variables u] and uf , 
that are uniformly distributed on ]0, 1[, by means 
of 



w, 



2 

- -log(w-) 



• Set 

ra(si. As) 



(As) 



(Al) 
(A2) 



(A3) 



^ sin[a(T^, + f)] 
[cos(K:)]^'/"^ 
cos [F, -a(y, + f)] 



In case of a = 1, Ti{si,As) — As is recovered. If 
adequate random number generators are applied for gen- 
eration of the variables Vi and Wi , the resulting random 
numbers are uncorrelated. From figure [7] it becomes 
evident, that the desired skewed a-stable Levy pdf ^ is 
matched. This algorithm therefore can be applied for ef- 
ficient numerical simulation of the process t{s) according 
to equation (|12b[) . 
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